投票率55.93%はいかがですか?

Primitive Analysis

作者
所属

Kansai University

公開

2023年2月9日

1 セットアップ

コード
pacman::p_load(tidyverse, readxl, naniar, marginaleffects, modelsummary)

raw_df <- read_xlsx("Data/Data.xlsx")

df <- raw_df %>%
  replace_with_na(replace = list(Q15.1 = c(88, 99),
                                 Q15.2 = c(88, 99),
                                 Q16.1 = c(88, 99),
                                 Q16.2 = c(88, 99),
                                 Q2.1  = 99,
                                 Q2.3  = 99,
                                 Q2.4  = 99,
                                 Q3_1  = 99,
                                 Q4.1  = 99,
                                 Q7.1  = 99,
                                 Q13.1 = c(6, 99))) %>%
  mutate(Q15.1 = 6 - Q15.1,
         Q15.2 = recode(Q15.2,
                        "1" = "5",
                        "2" = "4",
                        "3" = "2",
                        "4" = "1",
                        "5" = "3"),
         Q16.1 = 6 - Q16.1,
         Q16.2 = recode(Q16.2,
                        "1" = "5",
                        "2" = "4",
                        "3" = "2",
                        "4" = "1",
                        "5" = "3"),
         Female = sex - 1,
         Age    = age,
         Educ   = Q2.1,
         Income = Q2.3,
         Urban  = 5 - Q2.4,
         Effi   = Q3_1,
         Interest = 5 - Q4.1,
         Independent = if_else(Q7.1 == 11, 1, 0),
         Voted       = if_else(Q13.1 <= 3, 1, 0)) %>%
  rename(Treatment = votingrate) %>%
  filter(Q3_8 == 4)

df2 <- df %>%
  group_by(Treatment) %>%
  summarise(Outcome = mean(Q16.1, na.rm = TRUE),
            .groups = "drop")

2 55.93%に対する評価

コード
fit1 <- lm(Q15.1 ~ Female + Age + Educ + Income + Urban + 
             Effi + Interest + Independent + Voted,
   data = df)

modelsummary(fit1,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)
 (1)
(Intercept) 2.986 (0.174) ***
Female 0.141 (0.049) **
Age −0.003 (0.002)
Educ 0.000 (0.026)
Income 0.004 (0.005)
Urban 0.001 (0.026)
Effi −0.146 (0.022) ***
Interest −0.133 (0.038) ***
Independent −0.089 (0.049) +
Voted −0.247 (0.071) ***
Num.Obs. 1723
R2 0.093
R2 Adj. 0.088
AIC 4746.3
BIC 4806.2
Log.Lik. −2362.136
RMSE 0.95
コード
tidy(fit1, conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>%
  mutate(term = fct_inorder(term),
         term = fct_rev(term),
         sig  = if_else(p.value < 0.05, "sig", "insig")) %>%
  ggplot() +
  geom_vline(xintercept = 0) +
  geom_pointrange(aes(x = estimate, y = term,
                      xmin = conf.low, xmax = conf.high, color = sig)) +
  scale_color_manual(values = c("sig" = "black", "insig" = "gray70")) +
  labs(x = "係数(正の場合、ポジティブ評価)") +
  guides(color = "none") +
  theme_bw(base_size = 12) +
  theme(axis.title.y = element_blank())

3 実験(処置効果)

コード
fit2 <- lm(Q16.1 ~ poly(Treatment, 2) + 
             Female + Age + Educ + Income + Urban + 
             Effi + Interest + Independent + Voted, 
   data = df)

modelsummary(fit2,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)
 (1)
(Intercept) 2.399 (0.159) ***
poly(Treatment, 2)1 42.939 (1.023) ***
poly(Treatment, 2)2 16.364 (1.024) ***
Female 0.111 (0.045) *
Age 0.001 (0.001)
Educ 0.014 (0.025)
Income −0.003 (0.004)
Urban 0.006 (0.024)
Effi −0.104 (0.021) ***
Interest −0.034 (0.035)
Independent −0.021 (0.045)
Voted −0.015 (0.065)
Num.Obs. 1761
R2 0.544
R2 Adj. 0.541
AIC 4629.5
BIC 4700.6
Log.Lik. −2301.736
RMSE 0.89
コード
predictions(fit2, newdata = datagrid(Treatment = 20:80)) %>% 
  ggplot(aes(x = Treatment, y = estimate)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_point(data = df2, aes(x = Treatment, y = Outcome),
             color = "gray80", size = 3) +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率(%)", y = "投票率に対する評価 w/ 95% CI",
       caption = "注:点は投票率ごとの応答変数の平均値を示す。") +
  theme_bw(base_size = 12)

4 条件付き処置効果

4.1 性別

コード
fit3 <- lm(Q16.1 ~ poly(Treatment, 2) * Female + 
             Age + Educ + Income + Urban + 
             Effi + Interest + Independent + Voted, 
   data = df)

modelsummary(fit3,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)
 (1)
(Intercept) 2.396 (0.159) ***
poly(Treatment, 2)1 41.346 (1.439) ***
poly(Treatment, 2)2 17.152 (1.441) ***
Female 0.112 (0.045) *
Age 0.001 (0.001)
Educ 0.014 (0.025)
Income −0.003 (0.004)
Urban 0.007 (0.024)
Effi −0.104 (0.021) ***
Interest −0.032 (0.035)
Independent −0.022 (0.045)
Voted −0.018 (0.065)
poly(Treatment, 2)1 × Female 3.217 (2.049)
poly(Treatment, 2)2 × Female −1.580 (2.052)
Num.Obs. 1761
R2 0.545
R2 Adj. 0.541
AIC 4630.4
BIC 4712.5
Log.Lik. −2300.204
RMSE 0.89

予測値

コード
predictions(fit3,
            newdata = datagrid(Treatment = 20:80,
                               Female     = 0:1)) %>%
  mutate(Female = factor(Female, levels = 0:1, 
                         labels = c("男性", "女性"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Female), alpha = 0.25) +
  geom_line(aes(color = Female), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率(%)", y = "投票率に対する評価 w/ 95% CI", 
       color = "性別", fill = "性別") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果(処置変数)

コード
marginaleffects(fit3, variables = "Treatment",
                newdata = datagrid(Female = 0:1)) %>%
  mutate(Female = factor(Female, levels = 0:1, 
                         labels = c("男性", "女性"))) %>%
  ggplot(aes(x = Female, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  labs(x = "性別", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果(調整変数)

コード
marginaleffects(fit3, variables = "Female",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数(= 投票率)の値(%)", 
       y = "性別の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.2 年齢

コード
fit4 <- lm(Q16.1 ~ poly(Treatment, 2) * Age + 
             Female + Educ + Income + Urban + 
             Effi + Interest + Independent + Voted, 
   data = df)

modelsummary(fit4,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)
 (1)
(Intercept) 2.416 (0.158) ***
poly(Treatment, 2)1 31.986 (3.388) ***
poly(Treatment, 2)2 6.091 (3.351) +
Age 0.001 (0.001)
Female 0.110 (0.045) *
Educ 0.014 (0.024)
Income −0.003 (0.004)
Urban 0.006 (0.024)
Effi −0.107 (0.020) ***
Interest −0.035 (0.035)
Independent −0.022 (0.045)
Voted −0.024 (0.065)
poly(Treatment, 2)1 × Age 0.216 (0.063) ***
poly(Treatment, 2)2 × Age 0.204 (0.063) **
Num.Obs. 1761
R2 0.550
R2 Adj. 0.546
AIC 4611.3
BIC 4693.5
Log.Lik. −2290.673
RMSE 0.89

予測値

コード
predictions(fit4,
            newdata = datagrid(Treatment = 20:80,
                               Age       = seq(20, 80, by = 20))) %>%
  mutate(Age = factor(Age, levels = seq(20, 80, by = 20), 
                      labels = paste0(seq(20, 80, by = 20), "歳"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Age), alpha = 0.25) +
  geom_line(aes(color = Age), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率(%)", y = "投票率に対する評価 w/ 95% CI", 
       color = "年齢", fill = "年齢") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果(処置変数)

コード
marginaleffects(fit4, variables = "Treatment",
                newdata = datagrid(Age = 18:79)) %>%
  ggplot(aes(x = Age, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "年齢(歳)", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果(調整変数)

コード
marginaleffects(fit4, variables = "Age",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数(= 投票率)の値(%)", 
       y = "年齢の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.3 教育水準

コード
fit5 <- lm(Q16.1 ~ poly(Treatment, 2) * Educ + 
             Female + Age + Income + Urban + 
             Effi + Interest + Independent + Voted, 
   data = df)

modelsummary(fit5,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)
 (1)
(Intercept) 2.401 (0.159) ***
poly(Treatment, 2)1 42.410 (3.685) ***
poly(Treatment, 2)2 19.000 (3.745) ***
Educ 0.014 (0.025)
Female 0.110 (0.045) *
Age 0.001 (0.001)
Income −0.003 (0.004)
Urban 0.006 (0.024)
Effi −0.104 (0.021) ***
Interest −0.035 (0.035)
Independent −0.021 (0.045)
Voted −0.014 (0.065)
poly(Treatment, 2)1 × Educ 0.159 (1.069)
poly(Treatment, 2)2 × Educ −0.785 (1.073)
Num.Obs. 1761
R2 0.544
R2 Adj. 0.541
AIC 4632.9
BIC 4715.0
Log.Lik. −2301.453
RMSE 0.89

予測値

コード
predictions(fit5,
            newdata = datagrid(Treatment = 20:80,
                               Educ      = 1:5)) %>%
  mutate(Educ = factor(Educ, levels = 1:5, 
                       labels = c("中学校", "高校",
                                  "短大・高専・専門学校",
                                  "大学", "大学院"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Educ), alpha = 0.25) +
  geom_line(aes(color = Educ), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率(%)", y = "投票率に対する評価 w/ 95% CI", 
       color = "最終学歴(在学中を含む)", 
       fill = "最終学歴(在学中を含む)") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果(処置変数)

コード
marginaleffects(fit5, variables = "Treatment",
                newdata = datagrid(Educ = 1:5)) %>%
  ggplot(aes(x = Educ, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  scale_x_continuous(breaks = 1:5,
                     labels = c("中学校", "高校",
                                "短大・高専・専門学校",
                                "大学", "大学院")) +
  labs(x = "最終学歴(在学中を含む)", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果(調整変数)

コード
marginaleffects(fit5, variables = "Educ",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数(= 投票率)の値(%)", 
       y = "最終学歴の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.4 世帯収入

コード
fit6 <- lm(Q16.1 ~ poly(Treatment, 2) * Income + 
             Female + Age + Educ + Urban + 
             Effi + Interest + Independent + Voted, 
   data = df)

modelsummary(fit6,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)
 (1)
(Intercept) 2.399 (0.159) ***
poly(Treatment, 2)1 43.191 (1.615) ***
poly(Treatment, 2)2 17.810 (1.615) ***
Income −0.003 (0.004)
Female 0.112 (0.045) *
Age 0.001 (0.001)
Educ 0.015 (0.025)
Urban 0.006 (0.024)
Effi −0.104 (0.021) ***
Interest −0.034 (0.035)
Independent −0.022 (0.045)
Voted −0.012 (0.065)
poly(Treatment, 2)1 × Income −0.034 (0.193)
poly(Treatment, 2)2 × Income −0.222 (0.193)
Num.Obs. 1761
R2 0.544
R2 Adj. 0.541
AIC 4632.1
BIC 4714.2
Log.Lik. −2301.037
RMSE 0.89

予測値

コード
predictions(fit6,
            newdata = datagrid(Treatment = 20:80,
                               Income    = c(1, 10, 21))) %>%
  mutate(Income = factor(Income, levels = c(1, 10, 21))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Income), alpha = 0.25) +
  geom_line(aes(color = Income), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率(%)", y = "投票率に対する評価 w/ 95% CI", 
       color = "世帯収入(万円)", fill = "世帯収入(万円)") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果(処置変数)

コード
marginaleffects(fit6, variables = "Treatment",
                newdata = datagrid(Income = 1:21)) %>%
  ggplot(aes(x = Income, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "世帯収入", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果(調整変数)

コード
marginaleffects(fit6, variables = "Income",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数(= 投票率)の値(%)", 
       y = "世帯収入の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.5 居住地域

コード
fit7 <- lm(Q16.1 ~ poly(Treatment, 2) * Urban + 
             Female + Age + Educ + Income + 
             Effi + Interest + Independent + Voted, 
   data = df)

modelsummary(fit7,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)
 (1)
(Intercept) 2.398 (0.159) ***
poly(Treatment, 2)1 38.079 (3.507) ***
poly(Treatment, 2)2 14.729 (3.467) ***
Urban 0.007 (0.024)
Female 0.112 (0.045) *
Age 0.001 (0.001)
Educ 0.014 (0.025)
Income −0.003 (0.004)
Effi −0.104 (0.021) ***
Interest −0.035 (0.035)
Independent −0.023 (0.045)
Voted −0.015 (0.065)
poly(Treatment, 2)1 × Urban 1.622 (1.122)
poly(Treatment, 2)2 × Urban 0.539 (1.114)
Num.Obs. 1761
R2 0.545
R2 Adj. 0.541
AIC 4631.1
BIC 4713.3
Log.Lik. −2300.573
RMSE 0.89

予測値

コード
predictions(fit7,
            newdata = datagrid(Treatment = 20:80,
                               Urban     = 1:4)) %>%
  mutate(Urban = factor(Urban, levels = 1:4, 
                        labels = c("町村", "人口10万人未満の市",
                                   "人口10万人以上の市",
                                   "政令市・23区"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Urban), alpha = 0.25) +
  geom_line(aes(color = Urban), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率(%)", y = "投票率に対する評価 w/ 95% CI", 
       color = "居住地域", fill = "居住地域") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果(処置変数)

コード
marginaleffects(fit7, variables = "Treatment",
                newdata = datagrid(Urban = 1:4)) %>%
  ggplot(aes(x = Urban, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  scale_x_continuous(breaks = 1:4,
                     labels = c("町村", "人口10万人未満の市",
                                "人口10万人以上の市", "政令市・23区")) +
  labs(x = "最終学歴(在学中を含む)", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果(調整変数)

コード
marginaleffects(fit7, variables = "Urban",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数(= 投票率)の値(%)", 
       y = "居住地域の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.6 投票義務感

コード
fit8 <- lm(Q16.1 ~ poly(Treatment, 2) * Effi + 
             Female + Age + Educ + Income + Urban + 
             Interest + Independent + Voted, 
   data = df)

modelsummary(fit8,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)
 (1)
(Intercept) 2.402 (0.159) ***
poly(Treatment, 2)1 45.485 (3.337) ***
poly(Treatment, 2)2 12.301 (3.313) ***
Effi −0.104 (0.021) ***
Female 0.112 (0.045) *
Age 0.001 (0.001)
Educ 0.014 (0.025)
Income −0.003 (0.004)
Urban 0.006 (0.024)
Interest −0.033 (0.035)
Independent −0.022 (0.045)
Voted −0.018 (0.065)
poly(Treatment, 2)1 × Effi −0.672 (0.830)
poly(Treatment, 2)2 × Effi 1.070 (0.827)
Num.Obs. 1761
R2 0.545
R2 Adj. 0.541
AIC 4631.3
BIC 4713.4
Log.Lik. −2300.626
RMSE 0.89

予測値

コード
predictions(fit8,
            newdata = datagrid(Treatment = 20:80,
                               Effi      = c(1, 3, 5))) %>%
  mutate(Effi = factor(Effi, levels = c(1, 3, 5), 
                       labels = c("そう思う", "どちらともいえない",
                                  "まったくそう思わない"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Effi), alpha = 0.25) +
  geom_line(aes(color = Effi), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率(%)", y = "投票率に対する評価 w/ 95% CI", 
       color = "自分一人くらい投票しなくてもかまわない", 
       fill = "自分一人くらい投票しなくてもかまわない") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果(処置変数)

コード
marginaleffects(fit8, variables = "Treatment",
                newdata = datagrid(Effi = 1:5)) %>%
  mutate(Effi = factor(Effi, levels = 1:5, 
                       labels = c("そう思う", 
                                  "ややそう思う",
                                  "どちらともいえない",
                                  "あまりそう思わない",
                                  "まったくそう思わない"))) %>%
  ggplot(aes(x = Effi, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  labs(x = "選挙では大勢が投票するものだから、自分一人くらい投票しなくてもかまわない", 
       y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果(調整変数)

コード
marginaleffects(fit8, variables = "Effi",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数(= 投票率)の値(%)", 
       y = "投票義務感の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.7 政治関心

コード
fit9 <- lm(Q16.1 ~ poly(Treatment, 2) * Interest + 
             Female + Age + Educ + Income + Urban + 
             Effi + Independent + Voted, 
   data = df)

modelsummary(fit9,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)
 (1)
(Intercept) 2.405 (0.159) ***
poly(Treatment, 2)1 44.114 (4.195) ***
poly(Treatment, 2)2 6.940 (4.335)
Interest −0.031 (0.035)
Female 0.110 (0.045) *
Age 0.001 (0.001)
Educ 0.015 (0.025)
Income −0.003 (0.004)
Urban 0.005 (0.024)
Effi −0.104 (0.021) ***
Independent −0.018 (0.045)
Voted −0.021 (0.065)
poly(Treatment, 2)1 × Interest −0.391 (1.448)
poly(Treatment, 2)2 × Interest 3.331 (1.490) *
Num.Obs. 1761
R2 0.545
R2 Adj. 0.542
AIC 4628.4
BIC 4710.5
Log.Lik. −2299.194
RMSE 0.89

予測値

コード
predictions(fit9,
            newdata = datagrid(Treatment = 20:80,
                               Interest  = c(1, 4))) %>%
  mutate(Interest = factor(Interest, levels = c(1, 4), 
                           labels = c("まったく関心がない", 
                                      "とても関心がある"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Interest), alpha = 0.25) +
  geom_line(aes(color = Interest), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率(%)", y = "投票率に対する評価 w/ 95% CI", 
       color = "政治関心", fill = "政治関心") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果(処置変数)

コード
marginaleffects(fit9, variables = "Treatment",
                newdata = datagrid(Interest = 1:4)) %>%
  mutate(Interest = factor(Interest, levels = 1:4, 
                           labels = c("まったく関心がない", 
                                      "あまり関心がない",
                                      "ある程度関心がある",
                                      "とても関心がある"))) %>%
  ggplot(aes(x = Interest, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  labs(x = "政治関心", 
       y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果(調整変数)

コード
marginaleffects(fit9, variables = "Interest",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数(= 投票率)の値(%)", 
       y = "政治関心の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.8 支持政党

コード
fit10 <- lm(Q16.1 ~ poly(Treatment, 2) * Independent + 
             Female + Age + Educ + Income + Urban + 
             Effi + Interest + Voted, 
   data = df)

modelsummary(fit10,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)
 (1)
(Intercept) 2.396 (0.159) ***
poly(Treatment, 2)1 41.109 (1.334) ***
poly(Treatment, 2)2 16.945 (1.326) ***
Independent −0.020 (0.045)
Female 0.109 (0.045) *
Age 0.001 (0.001)
Educ 0.015 (0.025)
Income −0.003 (0.004)
Urban 0.005 (0.024)
Effi −0.104 (0.021) ***
Interest −0.035 (0.035)
Voted −0.012 (0.065)
poly(Treatment, 2)1 × Independent 4.424 (2.077) *
poly(Treatment, 2)2 × Independent −1.396 (2.085)
Num.Obs. 1761
R2 0.545
R2 Adj. 0.542
AIC 4628.5
BIC 4710.6
Log.Lik. −2299.230
RMSE 0.89

予測値

コード
predictions(fit10,
            newdata = datagrid(Treatment   = 20:80,
                               Independent = 0:1)) %>%
  mutate(Independent = factor(Independent, levels = 0:1, 
                              labels = c("あり", "なし"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Independent), alpha = 0.25) +
  geom_line(aes(color = Independent), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率(%)", y = "投票率に対する評価 w/ 95% CI", 
       color = "支持政党", fill = "支持政党") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果(処置変数)

コード
marginaleffects(fit10, variables = "Treatment",
                newdata = datagrid(Independent = 0:1)) %>%
  mutate(Independent = factor(Independent, levels = 0:1, 
                              labels = c("あり", "なし"))) %>%
  ggplot(aes(x = Independent, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  labs(x = "支持政党", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果(調整変数)

コード
marginaleffects(fit10, variables = "Independent",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数(= 投票率)の値(%)", 
       y = "支持政党の有無の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

4.9 投票有無

コード
fit11 <- lm(Q16.1 ~ poly(Treatment, 2) * Voted + 
             Female + Age + Educ + Income + Urban + 
             Effi + Interest + Independent, 
   data = df)

modelsummary(fit11,
             estimate  = c("{estimate} ({std.error}) {stars}"),
             statistic = NULL)
 (1)
(Intercept) 2.418 (0.159) ***
poly(Treatment, 2)1 45.016 (2.351) ***
poly(Treatment, 2)2 10.539 (2.486) ***
Voted −0.024 (0.065)
Female 0.110 (0.045) *
Age 0.001 (0.001)
Educ 0.014 (0.025)
Income −0.003 (0.004)
Urban 0.006 (0.024)
Effi −0.105 (0.021) ***
Interest −0.034 (0.035)
Independent −0.021 (0.045)
poly(Treatment, 2)1 × Voted −2.530 (2.615)
poly(Treatment, 2)2 × Voted 7.003 (2.729) *
Num.Obs. 1761
R2 0.546
R2 Adj. 0.543
AIC 4626.0
BIC 4708.1
Log.Lik. −2297.996
RMSE 0.89

予測値

コード
predictions(fit11,
            newdata = datagrid(Treatment = 20:80,
                               Voted     = 0:1)) %>%
  mutate(Voted = factor(Voted, levels = 0:1, 
                        labels = c("棄権", "投票"))) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_vline(xintercept = 55.93, linetype = 2, color = "red") +
  geom_ribbon(aes(fill = Voted), alpha = 0.25) +
  geom_line(aes(color = Voted), linewidth = 1) +
  coord_cartesian(ylim = c(1, 5)) +
  scale_x_continuous(breaks = c(20, 40, 55.93, 60, 80),
                     labels = c(20, 40, 55.93, 60, 80)) +
  scale_y_continuous(sec.axis = sec_axis(trans = ~.*1, name="",
                                         breaks = c(1:5),
                                         labels = c("低いと思う",
                                                    "やや低いと思う",
                                                    "ちょうどよいと思う",
                                                    "やや高いと思う",
                                                    "高いと思う"))) +
  labs(x = "刺激: 投票率(%)", y = "投票率に対する評価 w/ 95% CI", 
       color = "投票有無", fill = "投票有無") +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

限界効果(処置変数)

コード
marginaleffects(fit11, variables = "Treatment",
                newdata = datagrid(Voted = 0:1)) %>%
  mutate(Voted = factor(Voted, levels = 0:1, 
                        labels = c("棄権", "投票"))) %>%
  ggplot(aes(x = Voted, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_pointrange(size = 1, linewidth = 1) +
  labs(x = "投票有無", y = "処置変数の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

限界効果(調整変数)

コード
marginaleffects(fit11, variables = "Voted",
                newdata = datagrid(Treatment = 20:80)) %>%
  ggplot(aes(x = Treatment, y = estimate, 
             ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0) +
  geom_ribbon(fill = "gray30", alpha = 0.3) +
  geom_line(linewidth = 1) +
  labs(x = "処置変数(= 投票率)の値(%)", 
       y = "投票有無の限界効果 w/ 95% CI") +
  theme_bw(base_size = 12)

5 RDD

 40%が投票率評価の「高低」の分岐点?

コード
rdd::RDestimate(Q16.1 ~ Treatment, cutpoint = 40, data = df) %>%
  summary()

Call:
rdd::RDestimate(formula = Q16.1 ~ Treatment, data = df, cutpoint = 40)

Type:
sharp 

Estimates:
           Bandwidth  Observations  Estimate  Std. Error  z value  Pr(>|z|) 
LATE       3.192      265           0.4530    0.12171     3.722    0.0001977
Half-BW    1.596      118           0.3500    0.12496     2.801    0.0050972
Double-BW  6.384      487           0.3235    0.09812     3.297    0.0009785
              
LATE       ***
Half-BW    ** 
Double-BW  ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F-statistics:
           F      Num. DoF  Denom. DoF  p      
LATE       3.815  3         261         0.02115
Half-BW    2.536  2         115         0.16720
Double-BW  4.325  3         483         0.01010
コード
rdd::RDestimate(Q16.1 ~ Treatment, cutpoint = 40, data = df) %>%
  plot()

 ちなみにLATEが統計的に有意である閾値は40%と45%

コード
25:75 %>%
  enframe(name = "ID", value = "Cutpoint") %>%
  mutate(Fit = map(Cutpoint, ~rdd::RDestimate(Q16.1 ~ Treatment,
                                              cutpoint = .x,
                                              data = df)),
         Est = map_dbl(Fit, ~.x$est[1]),
         SE  = map_dbl(Fit, ~.x$se[1]),
         LL  = Est + qnorm(0.025) * SE,
         UL  = Est + qnorm(0.975) * SE,
         p   = map_dbl(Fit, ~.x$p[1]),
         Sig = if_else(p < 0.05, "Sig", "Insig")) %>%
  ggplot() +
  geom_hline(yintercept = 0) +
  geom_pointrange(aes(x = Cutpoint, y = Est, ymin = LL, ymax = UL,
                      color = Sig)) +
  labs(y = "LATE w/ 95% CI") +
  scale_color_manual(values = c("Sig" = "black", "Insig" = "gray80")) +
  guides(color = "none") +
  theme_bw()

 むろん、この40%になにかの根拠があるわけではない。むしろ、change-in-slopeモデルで分析した場合、50%という結果も。

コード
pacman::p_load(cpop)

cpop(df2$Outcome, df2$Treatment)

cpop analysis with n = 61 and penalty (beta)  = 8.221748

1  changepoint detected at x = 
 50
fitted values : 
  x0       y0 x1       y1   gradient  intercept       RSS
1 20 1.113460 50 1.424349 0.01036297  0.9062005 0.4664509
2 50 1.424349 80 4.110136 0.08952625 -3.0519633 0.6733183

overall RSS = 1.139769
cost = 85.89115
コード
cpop(df2$Outcome, df2$Treatment) %>%
  plot()

 閾値を50%にしたRDDの場合、LATEは統計的に有意ではない。

コード
rdd::RDestimate(Q16.1 ~ Treatment, cutpoint = 50, data = df) %>%
  summary()

Call:
rdd::RDestimate(formula = Q16.1 ~ Treatment, data = df, cutpoint = 50)

Type:
sharp 

Estimates:
           Bandwidth  Observations  Estimate  Std. Error  z value  Pr(>|z|)   
LATE       3.36       234           -0.3896   0.3493      -1.115   0.2647     
Half-BW    1.68        82           -0.3333   0.2348      -1.420   0.1556     
Double-BW  6.72       459           -0.2948   0.2000      -1.474   0.1404     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

F-statistics:
           F      Num. DoF  Denom. DoF  p      
LATE       2.118  3         230         0.19747
Half-BW    1.359  2          79         0.52593
Double-BW  4.270  3         455         0.01093
コード
rdd::RDestimate(Q16.1 ~ Treatment, cutpoint = 50, data = df) %>%
  plot()